/**
 * @file quic_svd.hpp
 * @author Mudit Raj Gupta
 *
 * Implementation of QUIC SVD class to perform Matrix Factorization
 * on the given matrix.
 *
 */

#ifndef __MLPACK_CORE_OPTIMIZERS_QUIC_SVD_QUIC_SVD_IMPL_HPP
#define __MLPACK_CORE_OPTIMIZERS_QUIC_SVD_QUIC_SVD_IMPL_HPP

namespace mlpack {
namespace quicsvd {

QuicSVD::QuicSVD(double epsilon, double delta, double maxIterations) : 
    epsilon(epsilon),
    delta(delta),
    maxIterations(maxIterations)
{
  //Nothing to do
}

QuicSVD::QuicSVD()
{
  this->epsilon = 0.9;
  this->delta = 0.1;
  this->maxIterations = 1000;
}

void QuicSVD::Apply(arma::mat& U, arma::mat& sigma, arma::mat& V, arma::mat A)
{
CosineTree Nroot,N;
CosineTreeBuilder build;
arma::mat V;
double mcSqErr = arma::norm(A,"fro");
std::map<double,CosineTree> Q;
Nroot = build.CTNode(A);
Q.insert(std::pair<double,CosineTree&>(0.0,Nroot));
CosineTree left, right, *p, child;
*p = NULL;
double errC;
size_t count=0;
size_t sl,sg;
sg = Nroot.NumPoints();
while(mcSqErr <= epsilon*norm(A,"fro") || count < maxIterations)
{
  left = p;
  right = p;
  N = Q.rbegin()->second;
  Q.erase(Q.rbegin());
  builder.CTNodeSplit(N,left,right);  
  //Remove N's contributed basis vector from V
  for(size_t i = 1; i < 3; i++)
  {
    child = root.Child(i);
    V = MGS(V);
  } 
  for(size_t i = 1; i < 3; i++)
  {
    child = root.Child(i);
    sl = child.NumPoints();
    errC = MCSqError(child.Data(),V,sl);   
    Q.insert(errC, child);
  }
  mcSqErr = MCSqError(A,V,sg);
  count ++;
}
}

void QuicSVD::RowSamplesLS(arma::mat A, arma::mat& S, arma::vec& ps,size_t n)
{
  arma::vec probability = arma::zeros<arma::vec>(A.n_rows,1);
  LSSampling(A, probability);
  arma::vec partialSortProb = probability;
  std::nth_element(partialSortProb.begin(), partialSortProb.begin() + n, partialSortProb.end());
  size_t count = 0;
  for(size_t i=0;i<probability.n_rows||count<n;i++)
  {
    if(probability(i)<=partialSortProb(n))
    {
      S.col(count) = A.col(i);
      ps(count) = probability(i);
      count ++; 
    }
  }
}

double QuicSVD::LowBound(double mu, double var, size_t n, double delta)
{

}

double QuicSVD::MCSqError(arma::mat A, arma::mat V, size_t n, double delta)
{
arma::mat S = arma::zeros<arma::mat>(A.n_rows,n);
arma::vec ps = arma::zeros<arma::vec>(A.n_rows,1);
S = RowSamplesLS(A,S,ps,n);
arma::vec wgtMagSq;
double mu,variance,magSqLB;
size_t count;
for (count=0;count<n;count++)
  wgtMagSq(count) = (1.0/ps(count))*arma::norm((S.col(count)*V),"fro");
mu = arma::mean(wgtMagSq);
variance = arma::var(wgtMagSq);
magSqLB = LowBound(mu,variance,n,delta);
return norm(A,"fro") - magSqLB;
}


void QuicSVD::ExtractSVD (arma::mat& U, arma::mat sigma, arma::mat& V, arma::mat A, arma::mat VCap)
{
  arma::mat tmp;
  tmp = arma::dot(A,VCap);
  tmp = arma::dot(strans(tmp),tmp);
  arma::mat UApprox,VApprox;
  arma::vec sApprox;
  arma::svd(UApprox,sApprox,VApprox,tmp); 
  V = VCap*VApprox;
  sigma = sqrt(arma::diagmat(sApprox));
  U = (A*VCap)*VApprox*arma::pinv(sigma);
}

arma::mat QuicSVD::MGS(arma::mat V)
{
  arma::mat VNorm= zeros<arma::mat>(V.n_rows,V.n_cols);
  for (size_t i=0;i<V.n_cols;i++)
  {
    VNorm.col(i) = V.col(i);
    for (size_t j=0;j<i;j++)
      VNorm.col(i) -= dot(VNorm.col(i), VNorm.col(j))*VNorm.col(j);    
  }
  for (size_t i=0;i<VNorm.n_cols;i++)
    VNorm.col(i) = (1.0/EuclideanNorm(VNorm.col(i)))*VNorm.col(i); 
  return VNorm;
}

void QuicSVD::LSSampling(arma::mat A, arma::vec& probability)
{
  Log::Info<<"LSSampling"<<std::endl;
  //Saving the frobenious norm of the matrix
  double normA = arma::norm(A,"fro");
  //Calculating probability of each point to be sampled
  for (size_t i=0;i<A.n_rows;i++)
    probability(i) = arma::norm(A.row(i),"fro")/normA;
}

double QuicSVD::EuclideanNorm(arma::vec A)
{
   Log::Info<<"EuclideanNorm"<<std::endl;
   //Calculating the Euclidean Norm
   return sqrt(arma::sum(arma::square(A)));
}

void QuicSVD::Epsilon(double e)
{
  epsilon = e;
}

double QuicSVD::Epsilon()
{
  return epsilon;
}

void QuicSVD::Delta(double d)
{
  delta = d;
}

double QuicSVD::Delta()
{
  return delta;
}

void QuicSVD::MaxIterations(size_t m)
{
  maxIterations = m;
}

size_t QuicSVD::MaxIterations()
{
  return maxIterations;
}

}; // namespace quicsvd
}; // namespace mlpack

#endif
